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Quadrature moments method for the simulation 
of turbulent reactive flows 


By Veiikatramanan Raman, Heinz Pitsch and Rodney 0. Fox f 


1. Motivation and Objectives 

In recent years, computational fluid mechanics has become one of the primary tools for 
design and optimization of chemical reactors. With stringent environmental constraints, 
close control of product selectivity and an estimate of by-products are essential in suc- 
cessfully operating chemical plants. The fast throughput and enhanced mixing conditions 
offered by turbulent flows are increasingly exploited in chemical reactors. Viable simula- 
tion methods for such flows should be able to model the complex interaction of reaction 
and turbulent flow. All known reaction models used in these simulations follow a segre- 
gated approach where different techniques are used to solve the momentum and scalar 
transport equations. It is assumed that reaction affects fluid flow only through the change 
in density. Usually a variable density flow solver is utilized that accepts the density field 
from the scalar handler to correct the flow field. The scalar transport scheme uses the flow 
properties to evaluate the local density and this iterative procedure is used to advance 
the solution in time. The Eulerian solution technique that is commonly used in solving 
scalar transport equations inherently does not contain information about sub-grid level 
processes. The adverse effect of neglecting sub-grid scalar fluctuations in cases where 
the scalar evolves through non-linear rate expressions is well known. In combustion pro- 
cesses which are characterized by fast chemistry, use of flamelet model (Pitsch & Steiner 
2000) or conditional moment closure (CMC) (Bilger 1993) based on a conserved scalar 
is known to be quite accurate in making qualitative as well as quantitative predictions. 
Such a method obviates the need for solving multiple scalar transport equations and 
is not restricted by the time-scale of individual reactions in the chemistry mechanism. 
The flamelet model, like all other reaction models, requires the specification of a scalar 
dissipation rate and assumes the shape of the PDF at the sub-grid level. The first-order 
CMC method models the conditional mean of the reacting scalars as a function of mix- 
ture fraction and local flow conditions. But the model ignores any fluctuations about the 
conditional mean and may not be viable in slow chemistry regimes. Both the CMC and 
flamelet model assume that all the reactive scalars can be parameterized by a single con- 
served scalar. On the other end of the spectrum, the transported-PDF scheme computes 
the sub-grid scalar-PDF in terms of a set of delta- functions. This method can be used 
in tandem with a flow solver like those based on the Reynolds-Averaged Navier Stokes 
(RANS) equation or Large-Eddy Simulation (LES) scheme to model reaction (Muradoglu 
et al. 1999; Haworth & El Tahry 1991). Such a formulation computes the profiles of all 
the species involved in the flow and leads to a closed form for the reaction source term. 
However, a mixing model is needed to describe the sub-grid mixing process and has been 
the focus of study in PDF methods (Subramaniam & Pope 1999; Valino 1995; Tsai & Fox 
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1998). The fact that no convincing mixing model exists is one of the important drawbacks 
of the method. 

The transported-PDF still is an attractive scheme in that the reaction source term 
appears closed. The inhomogeneous transport equation for the sub-grid PDF is multi- 
dimensional and cannot be solved using Eulerian grid techniques. Particle based Monte- 
Carlo schemes used with realistic chemistry show good agreement with experimental 
results (Xu & Pope 2000; Masri & Pope 1990; Roekaerts 1991; Raman et al. 2003a). 
However, the handling of a large number of Lagrangian particles along with detailed 
chemistry can be computationally expensive. When using LES as the flow solver, such 
schemes can be outright intractable in practical flow configurations. In this work, an 
equivalent Eulerian version of the transported-PDF method is formulated for LES of 
reactive flows. 


2. Filtered Quadrature Scheme for Variable Density Flows 


For use in LES of reactive flows, a filtered density function (FDF) is defined (Gao 
& O’Brien 1993) that prescribes the sub-filter joint composition density function of the 
scalars. 


F L (if>; x, t) 
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where Fl is the mass density function. This definition implies that in variable density 
flows, the FDF is the mass weighted and spatially filtered fine-grain density (Jaberi et al. 
1999). 
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where p is the filtered density. Using the scalar transport equation, the FDF is shown to 
evolve in multi-dimensional space as: 
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where u, is the Favre filtered velocity component, ip is the composition array and <p 
is the Favre filtered mean composition array. The molecular and turbulent diffusivi- 
ties are denoted by 7 and 7 1 respectively. The mixing frequency is modeled as fl m = 
Cq (7 + 74) / ({p) 1 Aq), where Aq is the filter width. This high- dimensional equation 
is traditionally solved using Monte-Carlo schemes (Colucci et al. 1998; Raman et al. 
2003&). As mentioned before, the particle based Monte-Carlo schemes are computation- 
ally expensive for large grids and hence become intractable even for simple chemistry. 
The Direct Quadrature Method of Moments (DQMOM) (Marchisio & Fox 2003; Wang 
& Fox 2003) is introduced here in the context of LES to provide an alternate tractable 
scheme. 

Figure 1 shows the fine-grained PDF composed of a set of delta-functions approximat- 
ing an arbitrary shaped PDF at a given point in the flow. Any particle based solution 
to the FDF equation will yield such an approximation. The FDF plot can be considered 
as a plot of normalized weights of particles in the computational cell. The accuracy of 
the approximation depends, among other factors, on the number of approximating delta 
functions. The N a delta-functions are characterized by their positions (4>ai) and their 
heights ( w a i ). The FDF is considered solved if for a given number of delta functions, the 
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Figure 1 . Approximation of a PDF using a finite number of delta functions. (Left) 
Transported PDF method and (Right) DQMOM method. 


values of <p a i and w a i are known. The DQMOM model formulates transport equations 
for these quantities based on pre-determined number of so-called environments, where 
each environment corresponds to a single delta-peak. The source terms involved in these 
equations ensure that the DQMOM equations match the moment equations of the scalar. 


In the multi-environment DQMOM model, the FDF is assumed to be of the form: 


N 

F L = P^ WiS (if - </>,;) . 

i = 1 


(2.4) 


The individual transport equations for the weights and locations of the delta-functions 
can be derived by substituting Eq. 2.4 into Eq. 2.3 (Fox 2003). For a single scalar case, 
this leads to 
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where Gj = u>i(p i and R(ip\ x, t) contains the mixing and reaction terms. 5 indicates the 
m-th derivative of the delta function. The properties of 5 function are used in rewriting 
the derivatives in terms of the filtered quantities (Pope 2000) . 
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and 




(2.7) 


Following Fox (2003), Eq. 2.5 can be rewritten in terms of the transport equations for 
Wi and Gi, using a t and bi to denote the respective equations. 
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where Cj = (r + r t )(V0j) 2 . In deriving these equations, the only assumption made so 
far is that the shape of the PDF is approximated by a finite-set of delta functions. For 
a 7V-environment model, 2 N source terms need to be specified. It should be noted that 
if the number of scalars is more than one, several cross moments can also be specified. 
For details on the feasibility and choice of moments refer to Fox (2003) and Marchisio 
& Fox (2003). Here, only the pure moments, namely the mean and variance of a scalar 
are used in fixing the source terms. Multiplying Eq. 2.8 by ip m and integrating over the 
composition space yields: 


N N 

, \ —rm —771—1 J 

(1- m) 2_^(j) <P A 

i= 1 i=l 


= m(m — 


N 

1)£ (p jWiCi + R 771 • 
2=1 


(2.9) 


Similarly the reaction and mixing terms can be integrated to yield: 
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Using Eqs 2.9 and 2.10 we can obtain the source terms a; and bi by solving a set of non- 
linear algberaic equations. In this work, N is set to two for all simulations. To further 
simplify the equations, we can arbitrarily set a.j to zero (Fox 2003). We then use m = 1,2 
to obtain the other source terms. Using the above moment equations, the non-linear 
system for bi reduces to: 
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From the above set of equations, the source terms are determined to be: 
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(2.12) 
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As discussed earlier, the number of environments can be considered to be equivalent to 
the particle number density used in a Lagrangian simulation. This implies that as the 
number of environments is increased, the accuracy of the scheme should increase (Wang 
& Fox 2003). The single scalar DQMOM can be extended to a multi- variate case using 
similar derivation technique (Fox 2003). In determining the source terms, one of the 
important criteria is that the pre-multiplier of the b vector in Eq. 2.11 is invertible. It 
should also be noted that cq need not be set to zero and can be determined through a 
set of extended non-linear equations (Marchisio & Fox 2003). 


3. Numerics 

The numerical implementation for reactive flow simulations consists of two parts - 
the flow solver and the scalar handler. Here, we use a LES flow solver along with three 
different implementations of the scalar transport equation. The first method solves for 
the mean mixture fraction and reaction progress variable along with the variance of 
the mixture fraction using Eulerian transport equations for these quantities. The second 
method uses a Lagrangian Monte-Carlo scheme to solve the FDF of the mixture-fraction 
and reaction progress variable. The third implementation uses the DQMOM based 2- 
environment model to obtain the first and second moments of the scalars. In all these 
cases, the density field is obtained from the scalar handler and fed to the flow solver which 
corrects the flow according to reaction. To first establish the viability of the DQMOM 
model, we only consider reactions with no heat release in constant density incompressible 
flows and thus do not use the feedback loop. The flow solver provides the one-way transfer 
of velocity and turbulence fields to the scalar handler. The following subsections briefly 
explain the numerical implementations of each of the components. 

3.1. LES solver 

The second order LES scheme uses an energy conserving formulation for the momentum 
equations (Pierce 2001). The eddy viscosity and diffusivity are computed using dynamic 
Smagorinsky model (Moin et al. 1991). The LES scheme also solves for the scalar trans- 
port equation using a semi-implicit scheme. In the current work, scalar equations are 
solved for the mixture fraction Z and the reaction progress variable Y. The numerical 
implementation uses an upwind based QUICK scheme that is designed to reduce numer- 
ical oscillations but also leads to some amount of numerical dissipation. 

For the sake of comparison, the mixture fraction variance transport equation is also 
simulated using production and dissipation terms consistent with the FDF formulation. 
Since sub-grid variance is a very small quantity in LES simulations, numerical dissipation 
can further increase the errors in the computation. To overcome this problem, a Z — Z 2 
system is solved rather than the variance transport equation. The scalar equations are 
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obtained by integrating the FDF equation over mixture fraction space after multiplying 
by Z and Z 2 . 

^ + VpUZ = v(r + r T )vz ( 3 . 1 ) 

^- + VpU^ = v(r + r T )vz2-f) ro (z2-zz), ( 3 . 2 ) 

with f l m defined as before providing the time scale for scalar dissipation. This formulation 
ensures that in the limit of zero dissipation, the sub-grid variance is identical to the 
analytical solution: 

Z 7 2 = Z2^Z = Z(l-Z). (3.3) 

The inlet and boundary conditions for Z 2 are identical to the filtered mixture frac- 
tion equation. This ensures that the sub-grid variance is zero at both the inlet and any 
boundaries in the domain. The chemical source term for the progress-variable equation 
is modeled based on the filtered means of the scalars and neglecting any sub-grid fluc- 
tuations. This “laminar” assumption leads to the following transport equation for the 
reactive scalar: 

^ + VpUF = V (T + T t ) VF + pS(Y) . (3.4) 

The LES solution of the mixture fraction mean and variance equation provides an inde- 
pendent way of checking the transported PDF and DQMOM methods. 

3.2. Transported PDF solver 

In the Lagrangian implementation, the FDF transport equation is solved using an en- 
semble of notional particles. The hybrid composition-PDF scheme used here uses the 
flow fields from the LES solver (that was described in the previous section) to evolve 
stochastic particles that are evenly distributed in the entire computational domain. The 
evolution equations for the particles are obtained from Eq. 2.3 by using techniques sim- 
ilar to the R.ANS based hybrid method (Pope 1985; Colucci et al. 1998). The SDE’s are 
solved for a system of particles which represents the FDF in composition space. The 
particles move in physical space according to: 

dx* = (u+ i v(r + r t )^ dt + ^ 2 ^ + _J T ^ dW, 

where x* is the instantaneous position of the particle and dW represents the 
diffusion term. The motion in composition space is due to mixing and reaction. 

d0* = [O m (0 - 0*) + S(0)] dt, (3.6) 

where 0* is the notional particle composition. Here the composition vector 0 = [Z,Y] 
evolves using reaction source terms 5(0) = [0, S(Z, Y)\. The use of the IEM model with 
the appropriate dissipation rate ensures that the variance dissipation is consistent with 
the LES solver. Since typical LES simulation use millions of computational cells, even a 
particle number density of 10-20 particles per cell will be computationally expensive. To 
ensure tractability, several novel computational algorithms are implemented. A pointer- 
based storage structure is used for handling the particles. Conventional algorithms use 
sorting procedures to keep the particles in order but are not feasible for such large 
numbers. Here an integer numbering procedure is used that uses pointer based linking of 
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particles to cells. Cells store only the index of the particles in their domain. Ownership 
of particles is transferred to the destination cell when particles move across faces. Such 
integer operations are inexpensive, even though they tend to loose some efficiency with 
iterations as the particles occupy discontinuous locations in memory. Particle motion is 
handled using an element-to-element tracking procedure (Subramaniam & Haworth 2000; 
Raman et al. 2003 a). This algorithm splits the motion into a sequence of sub-steps that 
move individual particles from face-to-face with intermediate re-interpolation of particle 
velocity. Reflective boundaries can then be handled directly. In all the simulations carried 
out here, a nominal number density of 30 particles per cell is used. Particle splitting and 
merging is utilized to minimize the fluctuations in the number density. The details of the 
implementation and some preliminary results can be found in Raman et al. (2003 b). 

3.3. DQMOM implementation 

The two-environment version of the model solves for the scalars uq (w 2 
G 12 , G 21 , G 2 2 - The working set of equations can be written as 

1 + Vpu'uji = v (r + r T ) Vuq. 

at. 

^^ + V]mGai = V(T + r r )VG ai + 5 ai . (3.8) 

It can be seen that the evolution equation for the environment (Eq. 3.7) is identical 
to the mixture-fraction equation with no source term. The weighted scalar equations 
(Eq. 3.8) contain source terms that can be decomposed into mixing, reaction and correc- 
tion components. These terms can be specified by first determining the location of the 
delta-functions in composition space in a given computational cell. 

0oi=— ■ (3-9) 

Wi 

The mixing source term takes the form: 

l = -S£i = Ki G a2 - w a2 G al ) . (3.10) 

The correction term is obtained using the location of the delta peaks as 

i 2 

s c a i = -SI 2 = = £(Pt + T) (V0 Qi ) 2 . (3.11) 

Yal ra2 j = i 

The reaction source term is computed based on the composition vector (/> i in each envi- 
ronment as 

S r ai = S a (&) . (3.12) 

The source term for the scalar equation Gai is then given by 

S ai = S™+S c ai + WiS r ai , (3.13) 

where the reaction source term is multiplied by the weight in order to be consistent with 
the formulation. 

The main implementation issue comes for the correction term. The correction term, 
compensates for the excess moment source term arising from the finite-peak repre- 
sentation of the FDF. As seen above, this source term involves the inverse of the separa- 
tion distance of the delta-peaks composition space. In regions of near-complete mixing, 


— 1 — uq), Gn, 
(3.7) 
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Figure 2. Instantaneous mixture fraction (left) and progress variable (right) contours 
simulated using the FDF scheme. 


this term will approach zero leading to very large correction terms. Several alternatives 
have been suggested (Wang & Fox 2003) but in this work, an ad-hoc limit is used such 
that for all points where the difference <p a i — (j > a 2 is less than e, this term is set to e. Such 
an implementation using e = 10 -4 was found to be stable for all the cases studied. 

Another important aspect to note is that the scalar solver does not maintain the 
bounds on the scalars. This might lead to peak locations outside the accessible range for 
the scalars. It was found that any attempt to limit these values led to a sharp decrease in 
the estimate of the variance. To counter this problem, for all grid cells that contain peaks 
outside the natural limits, the source terms were set to zero. The presence of reaction 
source terms with stiff kinetics was found to adversely affect this occurrence. Smaller 
time-steps that are determined based on reaction-time scale were found to alleviate this 
problem. 


4. Simulations 

The shear flow geometry of Mungal & Dimotakis (1984) is used to test the new scheme. 
The configuration consists of a planar shear layer formed by two streams entering at 
8.8 m/s and 22 m/s velocity. Though the experiment involves a low heat-release fast 
chemistry, in this work we have not implemented this mechanism. Since the purpose is to 
compare different reaction models, a simple first order mechanism of the type A + B —> P 
is used for modeling reaction. The rate expression is simplified using a mixture fraction- 
progress variable approach. Three different rate expressions that commonly occur in 
reacting flows are tested. The transported PDF scheme is also simulated for the same 
flow conditions and the results are compared with the DQMOM and LES simulations. 
A 256 x 128 grid spanning 80Z? in the axial direction and 40 D in the cross-stream 
direction is used. Here D is set such that the observation point in the experiment is 
around 501?. The inlet velocity profiles are assumed to be laminar with flat profiles. A 
development region corresponding to 10Z? extends into the domain where the two streams 
are separated by a splitter plate. 
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In all the simulations, the rate expression for the progress variable is of the form 


7'y = K 



1 - Z 
1 - Z Kt 



(4.1) 


The rate constant K is varied to study the effect of non-linearity on DQMOM predic- 
tions. In the first case, the rate constant is set to a value of I\ = 2. Figure 2 shows 
the instantaneous mixture fraction and progress variable contours near the center of the 
domain. It shows the vortex like structure common to shear flows and the presence of 
highly mixed reactants at the center of such vortices. The peak in source term and mean 
progress variable are observed near the centerline. The vortices were found to stretch 
to a maximum of 10D in the cross-stream direction. It was found that the DQMOM 
model exhibits similar high reaction-rate zones. However, the LES simulations show a 
thick reaction zone with maximum allowable reaction rate at each location. In order to 
compare the steady-state trends, mean and variance of all scalars were time-averaged 
for at least one flow through time. Figure 3 shows the cross-stream profiles of the mix- 
ture fraction and sub-grid variance computed from all three schemes. Theoretically, the 
sub-grid variance obtained from all these methods should be identical. However, the dif- 
ferences in the implementation cause some discrepancy. Nevertheless, the time-averaged 
filtered mixture fraction and sub-grid variance predicted by all the schemes are in good 
agreement, thereby validating both the DQMOM and transported-PDF implementation. 


The time-averaged mean and sub-grid variance of the reaction progress-variable ob- 
tained from the different schemes show some interesting features (Fig. 4). The sub-grid 
variance is non-zero only for the DQMOM and transported-PDF schemes and is set to 
zero for the LES scheme. In this context the LES solver can be considered as a one- 
environment model with complete sub-filter mixing. If the transported-PDF scheme is 
considered as a multi-environment DQMOM model, the particle scheme with a nominal 
number density of N represents an N-environment decomposition of the FDF. The cross 
stream profiles of the mean show that the DQMOM method provides a vast improvement 
over the one-environment solution. The mean profile shows that in spite of the simple 
rate expression, second moment terms cannot be neglected. The differences between the 
LES and DQMOM models are highest in the initial section where the effect of unmixed 
reactants will be very important. Since the inflow is laminar, the mixing layer itself does 
not become turbulent until about X = 7.5. However, the LES solver predicts very high 
reaction rates in even these laminar regions where low mixing should essentially keep 
the reaction rates to a very low value. This “early-ignition” is observed in the profiles at 
X = 20 where the mean value predicted by the LES solver is atleast 50% higher than 
that predicted by the DQMOM model. Surprisingly, the sub-grid variance profile from 
the DQMOM scheme also shows very good agreement with the PDF scheme. This essen- 
tially implies that the third and higher moments of the reactive scalar can be neglected 
for this chemistry scheme. 

In the next case, a more complex rate expression is implemented. Reaction rates ap- 
pearing in combustion have a strong dependence on temperature. Most source terms 
have an exponential dependence on local temperature. To simulate such a condition, the 
reaction rate constant was set to 

if= 1000exp(— 1 _ > < 1 8 - i y ( y) ), (4.2) 

where b denotes the degree of dependence on temperature. For practical combustion 
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Figure 3. Comparison of time-averaged mean and variance of mixture fraction using a 
constant at different axial locations. ( ) DQMOM, ( ) LES and ( ) PDF. 

applications, b is usually set to values between 5 and 6. However, such high values lead 
to extinction for the present flow configuration. Instead, a lower value of 1 is chosen so 
that the reaction zone can be anchored near the splitter plate. This implies a weaker 
dependence on temperature but nevertheless makes the rate expression non-linear. Fig. 5 
shows the cross-stream profiles of reaction progress variable. Here again, it can be seen 
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Figure 4. Comparison of time-averaged (left) mean and (right) variance of reaction progress 
variable using a constant rate constant at (top) X=20, (middle) X=30 and (bottom) X=50. 
( ) DQMOM, ( ) LES and ( ) PDF. 

that the DQMOM predictions of both the mean and the variance of the reactive scalar 
are in good agreement with the transported-PDF results. The LES predictions show fast 
rates consistent with the laminar assumption. It should be noted that depending on the 
reaction rates, the complete mixing assumption can also lead to quenching. 
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The third and final case uses a polynomial rate function which is common to chemical 
engineering assumptions. Reduced chemical mechanisms like the chloromethane reactions 
(West et al. 1999) may even involve non-integer moments of the scalar variable. Here, 
the rate is defined as 


r Y 


= 1000 


0.0001 + F 2 
1.0 + Y 


(4.3) 


This expression ensures that the reaction proceeds without the need for an ignition source. 
The predicted mean and variance (Fig. 6) indicate good agreement with the transported- 
PDF scheme. This clearly shows that the DQMOM model with just two environments 
is able to drastically improve the single-environment predictions obtained from the LES 
solver. 

In terms of computational requirements, the particle based transported-PDF solver 
was nearly 5 times slower than the DQMOM scheme even for such simple geometries. 
For complex configurations, the memory requirements of a large ensemble of particles 
can further slow down the simulation. However, the DQMOM scheme it not without 
limitations. In particular, the time-scale limitations of a stiff-chemistry source term can 
reduce the time-step used in the DQMOM model while the robust chemistry solvers 
(e.g ISAT (Pope 1997)) can increase the computational speed of the transported-PDF 
model. In addition, detailed chemistry mechanisms will require a large number of scalar 
transport equations and may eventually diminish the advantages over the particle scheme. 
In spite of these observations, present study shows that the DQMOM scheme is a viable 
altrernative to the Monte-Carlo based transported-PDF model and needs to be further 
explored for multi-species systems. 


5. Conclusions and Future Work 

A sub-filter model for reactive flows, namely the DQMOM model, was formulated for 
LES using the filtered mass density function. Transport equations required to determine 
the location and size of the delta-peaks were then formulated for a 2-peak decomposition 
of the FDF. The DQMOM scheme was implemented in an existing structured-grid LES 
solver. Simulations of scalar shear layer using an experimental configuration showed that 
the first and second moments of both reactive and inert scalars are in good agreement 
with a conventional Lagrangian scheme that evolves the same FDF. Comparisons with 
LES simulations performed using laminar chemistry assumption for the reactive scalar 
show that the new method provides vast improvements at minimal computational cost. 

Currently, the DQMOM model is being implemented for use with the progress vari- 
able/mixture fraction model of Pierce (2001). Comparisons with experimental results and 
LES simulations using a single-environment for the progress- variable are planned. Future 
studies will aim at understanding the effect of increase in environments on predictions. 
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Figure 6. Comparison of time-averaged (left) mean and (right) variance of reaction progress 

variable for polynomial reaction rate at (top) X=30 and (bottom) X=50. ( ) DQMOM, 

( ) LES and ( ) PDF. 
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